In [1]:
from query import *
import panel as pn
import geoviews as gv
import geoviews.feature as gf
import plotly.express as px
import geopandas as gpd
from shapely.geometry import Polygon
gv.extension("bokeh")
VARIABLE = "2m_temperature"
MIN_LAT = 55
MAX_LAT = 85
MIN_LON = -75
MAX_LON = -10
START_DATETIME = "2023-01-01 00:00:00"
END_DATETIME = "2023-01-03 09:00:00"
1. Single Value Aggregation¶
Get a single aggregation min/max/average value in the certain spatial area during a certain period of time
E.g., get the average temperature in the box area during the period of time.
In [2]:
single_value_aggregation_query(
variable=VARIABLE,
min_lat=MIN_LAT,
max_lat=MAX_LAT,
min_lon=MIN_LON,
max_lon=MAX_LON,
start_datetime=START_DATETIME,
end_datetime=END_DATETIME,
aggregation_method="mean"
)
Out[2]:
<xarray.Dataset> Size: 4B
Dimensions: ()
Data variables:
t2m float32 4B 256.22. Time series query¶
Get the time series for a variable's aggregation value in a certain spatial area during a certain period of time
E.g., get the daily average temperature in the box area during the period of time.
In [3]:
ts = time_series_query(
variable=VARIABLE,
min_lat=MIN_LAT,
max_lat=MAX_LAT,
min_lon=MIN_LON,
max_lon=MAX_LON,
start_datetime=START_DATETIME,
end_datetime=END_DATETIME,
time_resolution="day",
aggregation_method="mean"
)
ts
Out[3]:
<xarray.Dataset> Size: 36B
Dimensions: (time: 3)
Coordinates:
* time (time) datetime64[ns] 24B 2023-01-01 2023-01-02 2023-01-03
Data variables:
t2m (time) float32 12B 255.3 256.7 257.22.1 visualize time series¶
In [4]:
ts['t2m'].plot()
Out[4]:
[<matplotlib.lines.Line2D at 0x12f5c3c50>]
3. Heatmap query (single layer)¶
Get a single-layer, aggregated heatmap of a variable at a certain spatial area in a certain period of time
E.g., get a heatmap for max temperature for the box area during the period of time
In [5]:
hm = heat_map_query_single_layer(
variable=VARIABLE,
min_lat=MIN_LAT,
max_lat=MAX_LAT,
min_lon=MIN_LON,
max_lon=MAX_LON,
start_datetime=START_DATETIME,
end_datetime=END_DATETIME,
spatial_resolution=0.5,
spatial_agg_method="max",
)
hm
Out[5]:
<xarray.Dataset> Size: 32kB
Dimensions: (latitude: 60, longitude: 130)
Coordinates:
* longitude (longitude) float32 520B -74.88 -74.38 -73.88 ... -10.88 -10.38
* latitude (latitude) float32 240B 84.88 84.38 83.88 ... 56.38 55.88 55.38
Data variables:
t2m (latitude, longitude) float32 31kB 250.2 250.3 ... 281.1 281.13.1 visualize single-layer heatmap¶
In [6]:
hm['t2m'].plot(cmap="coolwarm")
Out[6]:
<matplotlib.collections.QuadMesh at 0x12f62ffd0>
4. Heatmap query (multi layer)¶
Get a multi-layer heatmap of a variable at a certain spatial area in a certain period of time
Each layer represents a time during the period of time.
AKA. get a raster cube of a variable at a certain spatial area and spatial resolution, in a certain period of time and temporal resolution.
E.g., get raster cude of temperature in the box area (at 2degree * 2degree resolution) during the period of time (at daily resolution)
In [7]:
hmm = heat_map_query_multi_layer(
variable=VARIABLE,
min_lat=MIN_LAT,
max_lat=MAX_LAT,
min_lon=MIN_LON,
max_lon=MAX_LON,
start_datetime=START_DATETIME,
end_datetime=END_DATETIME,
spatial_resolution=2,
spatial_agg_method="mean",
time_resolution="day",
time_agg_method="mean",
)
hmm
Out[7]:
<xarray.Dataset> Size: 6kB
Dimensions: (time: 3, latitude: 15, longitude: 32)
Coordinates:
* longitude (longitude) float32 128B -74.12 -72.12 -70.12 ... -14.12 -12.12
* latitude (latitude) float32 60B 84.12 82.12 80.12 ... 60.12 58.12 56.12
* time (time) datetime64[ns] 24B 2023-01-01 2023-01-02 2023-01-03
Data variables:
t2m (time, latitude, longitude) float32 6kB 249.4 249.4 ... 282.4
Attributes:
Conventions: CF-1.6
history: 2024-04-10 18:52:20 GMT by grib_to_netcdf-2.25.1: /opt/ecmw...4.1 visualize multi-layer heatmap¶
In [8]:
gv_dataset = gv.Dataset(hmm)
gv_dataset.to(gv.Image)
gv_image = gv_dataset.to(gv.Image)
gv_image.opts(cmap="coolwarm", colorbar=True, tools=["hover"])
gv_image = gv_image * gf.coastline
gv_image.opts(width=700, height=500)
# pn.panel(gv_image, width=700, height=500).servable()
Out[8]:
5. Value criteria query¶
Get a raster cube where each cell satisfies the value criteria
E.g., get raster cube in the box area (at 2degree * 2degree resolution) during the period of time (at daily resolution), where the daily max temperature is > 255
In [9]:
vc = value_criteria_query(
variable=VARIABLE,
min_lat=MIN_LAT,
max_lat=MAX_LAT,
min_lon=MIN_LON,
max_lon=MAX_LON,
start_datetime=START_DATETIME,
end_datetime=END_DATETIME,
spatial_resolution=2,
spatial_agg_method="max",
time_resolution="day",
time_agg_method="max",
criteria_predicate="<",
criteria_value=255
)
vc
Out[9]:
<xarray.Dataset> Size: 6kB
Dimensions: (time: 3, latitude: 15, longitude: 32)
Coordinates:
* longitude (longitude) float32 128B -74.12 -72.12 -70.12 ... -14.12 -12.12
* latitude (latitude) float32 60B 84.12 82.12 80.12 ... 60.12 58.12 56.12
* time (time) datetime64[ns] 24B 2023-01-01 2023-01-02 2023-01-03
Data variables:
t2m (time, latitude, longitude) float32 6kB 251.2 251.3 ... nan nan
Attributes:
Conventions: CF-1.6
history: 2024-04-10 18:52:20 GMT by grib_to_netcdf-2.25.1: /opt/ecmw...5.1 visualize value criteria result¶
In [10]:
gv_dataset = gv.Dataset(vc)
gv_dataset.to(gv.Image)
gv_image = gv_dataset.to(gv.Image)
gv_image.opts(cmap="coolwarm", colorbar=True, tools=["hover"])
gv_image = gv_image * gf.coastline
gv_image.opts(width=700, height=500)
# pn.panel(gv_image, width=700, height=500).servable()
Out[10]:
6. Area finding query¶
Find spatial area that satisfies a certain criteria.
E.g., find area in the box area where every daily max temperature is < 255
In [11]:
af = area_finding_query(
variable=VARIABLE,
min_lat=MIN_LAT,
max_lat=MAX_LAT,
min_lon=MIN_LON,
max_lon=MAX_LON,
start_datetime=START_DATETIME,
end_datetime=END_DATETIME,
spatial_resolution=2,
spatial_agg_method="max",
time_resolution="day",
time_agg_method="max",
criteria_predicate="<",
criteria_value=255,
any_or_all="all"
)
af
Out[11]:
<xarray.Dataset> Size: 6kB
Dimensions: (time: 3, latitude: 15, longitude: 32)
Coordinates:
* longitude (longitude) float32 128B -74.12 -72.12 ... -14.12 -12.12
* latitude (latitude) float32 60B 84.12 82.12 80.12 ... 60.12 58.12 56.12
* time (time) datetime64[ns] 24B 2023-01-01 2023-01-02 2023-01-03
Data variables:
t2m (time, latitude, longitude) float32 6kB 251.2 251.3 ... nan
spatial_mask (latitude, longitude) bool 480B True True ... False False
Attributes:
Conventions: CF-1.6
history: 2024-04-10 18:52:20 GMT by grib_to_netcdf-2.25.1: /opt/ecmw...6.1 visualize area finding¶
In [12]:
df = af["spatial_mask"].to_dataframe().reset_index()
df["latitude"] = df["latitude"] - 1
df["longitude"] = df["longitude"] - 1
df["latitude2"] = df["latitude"] + 2
df["longitude2"] = df["longitude"] + 2
gdf = gpd.GeoDataFrame(df, geometry=[Polygon([(x, y), (x, y2), (x2, y2), (x2, y)]) for x, y, x2, y2 in zip(df["longitude"], df["latitude"], df["longitude2"], df["latitude2"])])
fig = px.choropleth_mapbox(
gdf,
geojson=gdf.geometry,
locations=gdf.index,
color="spatial_mask",
mapbox_style="carto-positron",
center={"lat": gdf["latitude"].mean(), "lon": gdf["longitude"].mean()},
zoom=1,
)
fig.show()
7. Time period query¶
Find time period that satisfies a certain criteria
E.g., find days during the time period in the box area where all daily max temperature is < 282
In [13]:
tp = time_period_query(
variable=VARIABLE,
min_lat=MIN_LAT,
max_lat=MAX_LAT,
min_lon=MIN_LON,
max_lon=MAX_LON,
start_datetime=START_DATETIME,
end_datetime=END_DATETIME,
spatial_resolution=2,
spatial_agg_method="max",
time_resolution="day",
time_agg_method="max",
criteria_predicate="<",
criteria_value=282,
any_or_all="all"
)
tp
Out[13]:
<xarray.Dataset> Size: 6kB
Dimensions: (time: 3, latitude: 15, longitude: 32)
Coordinates:
* longitude (longitude) float32 128B -74.12 -72.12 -70.12 ... -14.12 -12.12
* latitude (latitude) float32 60B 84.12 82.12 80.12 ... 60.12 58.12 56.12
* time (time) datetime64[ns] 24B 2023-01-01 2023-01-02 2023-01-03
Data variables:
t2m (time, latitude, longitude) float32 6kB 251.2 251.3 ... nan nan
time_mask (time) bool 3B True False False
Attributes:
Conventions: CF-1.6
history: 2024-04-10 18:52:20 GMT by grib_to_netcdf-2.25.1: /opt/ecmw...7.1 visualize time period finding¶
In [14]:
df = tp["time_mask"].to_dataframe().reset_index()
fig = px.line(df, x="time", y="time_mask")
fig.update_yaxes(autorange="reversed")
fig.show()
8. Shape query¶
In [15]:
shape_file_path = "../data/geoshape/greenland.geojson"
gdf_shape = gpd.read_file(shape_file_path)
gdf_shape.plot()
Out[15]:
<Axes: >
In [16]:
sq = shape_query(
variable=VARIABLE,
shape_fpath=shape_file_path,
start_datetime=START_DATETIME,
end_datetime=END_DATETIME,
)
sq
Out[16]:
<xarray.Dataset> Size: 5MB
Dimensions: (time: 58, latitude: 95, longitude: 244)
Coordinates:
* longitude (longitude) float32 976B -73.0 -72.75 -72.5 ... -12.5 -12.25
* latitude (latitude) float32 380B 83.5 83.25 83.0 ... 60.5 60.25 60.0
* time (time) datetime64[ns] 464B 2023-01-01 ... 2023-01-03T09:00:00
Data variables:
t2m (time, latitude, longitude) float32 5MB 248.7 248.8 ... 281.8
shape_mask (latitude, longitude) bool 23kB False False ... False False
Attributes:
Conventions: CF-1.6
history: 2024-04-10 18:52:20 GMT by grib_to_netcdf-2.25.1: /opt/ecmw...8.1 visualize shape query¶
In [17]:
# mask data
masked = sq["t2m"].where(sq["shape_mask"])
In [18]:
# visualize the mask
sq["shape_mask"].plot()
Out[18]:
<matplotlib.collections.QuadMesh at 0x13cc2ab10>
In [19]:
# visualize a layer
masked[0].plot()
Out[19]:
<matplotlib.collections.QuadMesh at 0x133c8ffd0>
In [20]:
# visualize all layers
gv_dataset = gv.Dataset(masked)
gv_dataset.to(gv.Image)
gv_image = gv_dataset.to(gv.Image)
gv_image.opts(cmap="coolwarm", colorbar=True, tools=["hover"])
gv_image = gv_image * gf.coastline
gv_image.opts(width=700, height=500)
Out[20]: